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1.  INTRODUCTION 

In  this  report  we  discuss  three  problems  associated  with  the 
numerical  solution  of  the  transonic  unsteady  small-dlsturbance  equation 

(1.1) 

For  transonic  aeroelastlc  and  flutter  calculations,  Weatherlll 
et  al.  (1975)  and  Traci  et  al.  (1975)  have  reported  a numerical  method 
for  computing  small  transonic  unsteady  harmonic  perturbations  to  a 
steady  flow.  The  potential  function  Is  expressed  In  a series  of  Increasing 
powers  of  a small  parameter  which  measures  the  amplitude  of  an  unsteady 
disturbance  to  the  boundary  (e.g.,  a thin  airfoil  undergoing  harmonic 
oscillation) . This  results  In  a sequence  of  boundary-value  problems  for 
the  perturbation  potentials.  The  zeroth  order  Is  just  the  steady 
problem  and  Is  solved  by  a type-dependent  finite-difference  scheme  with  | 

a line  relaxation  procedure.  The  first-order  problem  results  In  a j 

linear  equation  of  mixed  type  for  the  perturbation  potential.  It  Is  ; ■ 

i , 

solved  by  the  same  procedure.  j 

The  method  of  Weatherlll  and  Traci  Is  successful  for  only  j 

small,  reduced  frequencies.  Beyond  a critical  frequency  (for  a given  \ 

Mach  number),  the  relaxation  solutions  diverge.  | 

In  section  2 of  this  report  we  discuss  several  remedies  for  | : 

this  problem  and  present  some  computational  examples . j 

The  formulation  of  Weatherlll  and  Traci  does  not  account  for 
the  unsteady  motion  of  the  shock  wave.  In  section  3 of  this  report,  we 
discuss  a method  to  account  for  this,  and  we  present  a calculated  example. 

The  third  problem  we  have  studied  Is  the  solution  of  equation  i 

(1.1)  by  finite-element  methods.  An  extensive  review  of  finite-element 
and  finite-difference  methods  for  transonic  flow  was  reported  by  Hafez,  ! 

Wellford  and  Murman  (1977).*  A finite-element  method  for  solving  equation  | 

(1.1)  was  formulated,  and  It  Is  presented  In  section  4.  The  proposed 
method  may  be  used  for  both  steady  and  unsteady  flow  problems.  For  steady 


This  study  was  partially  sponsored  by  NASA-Langley  Research  Center  under 
contract  NASl-14246. 
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[ flow  problems,  the  time  t may  be  taken  as  an  artificial  time.  Calcu- 

j lations  for  the  steady  problem  are  given  in  section  4.  Calculations  for 

the  unsteady  problem  have  not  yet  been  completed. 
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2.  NUMERICAL  STABILITY  OF  THE  HARMONIC  APPROACH 

2. 1 Convergence  Analysis 

Weatherill  et  al.  (1975)  analyzed  the  convergence  of  the 
relaxation  procedure  for  the  perturbation  potential.*  The  finite- 
difference  approximation  results  In  a matrix  that  Is  not  always  positive 
definite  (depending  on  the  reduced  frequency).  Remedies  for  this  problem 
are  presented  In  this  section,  and  simple  Illustrative  examples  are 
computed.  Alternative  methods  for  solving  the  problem  are  suggested. 

The  Helmholtz  equation  was  suggested  by  Weatherill  et  al.  (1975) 
for  studying  the  divergence  of  relaxation  solutions  that  he  and  his 
colleagues  encountered  when  solving  the  unsteady  transonic  flow  equation. 

The  Helmholtz  equation 

+ 4i  + u)^(|>  >■  0 , 0 < X < L, 

XX  yy  1 

0 < y < L2  (2.1) 

with  (|)  given  on  the  boundaries  of  a rectangle  of  length  L2  and 
width  L^  , Is  simpler  than  the  original  equation,  but  produces  similar 
behavior  when  solving  It  by  relaxation  Is  attempted.  That  Is,  for 
values  of  OJ  greater  than  a specific  value  , line  relaxation  diverges. 

The  one-dlmenslonal  analog  to  Helmholtz's  equation  produces  the  same 
behavior,  yet  It  Is  easier  to  handle.  For  numerical  Illustration  we 
shiill  therefore  consider  the  boundary-value  problem 


This  linear  equation  Is  singular  at  the  sonic  line.  The  fundamental 
solutions  (and  the  Green's  functions)  are  different  In  the  subsonic, 
sonic,  and  supersonic  regions,  and  constructing  the  solution  along  these 
lines  Is  a formidable  task.  Use  of  a Fourier  transform  In  the  latter 
direction  to  reduce  the  two-dimensional  problem  to  a one-dlmenslonal  one 
will  Introduce  convolution  Integrals  (since  the  coefficients  K - 4)° 
and  are  functions  of  x and  y ).  Hence,  we  are  left  with  finite 

differences. 
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0 < X < 1 


(2.2a) 


(p  + u)^(p  ■ 0 

<t«(0)  - (2.2b) 

<«>(1)  - • (2- 2c) 

Defining  a computational  net  with  net  spaclngs  Ax  and  Ay 
on  the  rectangle 


0 < X < 

0 < y < L2 

and  approximating  the  derivatives  of  equation  (2.1)  by  central  difference 
formulas  leads  to  a system  of  algebraic  equations  of  the  form 

- f , (2.3) 

where  A Is  the  matrix  of  coefficients,  (|)  Is  a vector  for  the  unknowns 
at  the  grid  points,  and  f Is  a vector  for  the  nonhomogeneous  terms  and 
the  boundary  conditions. 

Trying  to  solve  equation  (2.3)  by  a line  relaxation  procedure, 
Weatherlll  et  al.  (1975)  observed  the  frequency-dependent  limitation  on 
the  convergence  of  the  relaxation  method.  They  showed  that  relaxation 
procedures  converge  If  and  only  If  A Is  positive  definite.  For  fixed 
values  of  , Ax  , and  Ay  , the  eigenvalues  of  A are  functions 

of  u , and  a critical  frequency  for  which  the  smallest  eigenvalue 

Is  zero  exists.  Above  this  value,  relaxation  procedures  diverge.  We 
describe  below  two  solutions  to  this  problem. 

2.2  Remedies  to  the  Divergence  Problem 

2.2.1  The  A*A  Method 

To  solve  eq.  (2.3)  by  a relaxation  procedure  that  converges 
regardless  of  the  value  of  u , we  must  modify  the  equation  so  that  the 
new  system  matrix  Is  positive  definite.  Thus,  we  multiply  both  sides  of 
eq.  (2.3)  by  A*  , the  conjugate  transpose  of  A , to  give 
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A*A(|)  « A*f 

(2.4) 

or 

where 

= A*A 

and 

- A*f  . 

is  positive  definite  except  when  O)  Is  the  natural  frequency.  In 
which  case  A Is  singular.  If  A Is  not  singular,  eqs.  (2.3)  and  (2.4) 
are  equivalent  and  have  the  same  solution.  Relaxation  procedures  will 
always  converge  for  the  new  system  eq.  (2.4)  (as  long  as  the  relaxation 
factor  Is  positive  and  less  than  2) . The  disadvantage  of  this  remedy  Is 
that  the  bandwidth  of  A*A  Is  almost  twice  that  of  A , so  that  the 
rate  of  convergence  may  be  slow.  In  particular.  If  A Is  Ill-conditioned, 
A*A  will  be  even  more  so.  For  high  frequencies  the  grid  size  must  be 
quite  small;  hence,  the  number  of  equations  Is  large,  and  the  eigenvalues 
are  close  to  each  other.  The  smallest  eigenvalue  of  A*A  will  be  very 
close  to  zero.  This  causes  the  rate  of  convergence  to  be  slow.* 

We  demonstrate  that  the  method  described  above  enables  us  to 
solve  the  Helmholtz  equation  for  values  of  u greater  than  by 

considering  the  one-dlmenslonal  analog  to  the  Helmholtz  equation  (2.2). 

To  solve  this  problem  numerically,  we  choose  a set  of  J discrete  points 
Xj  with  uniform  spacing  Ax  In  the  Interval  of  Interest  (0  < x < 1). 

We  write  the  finite-difference  approximation  to  equation  (2.2a)  at  each  of 
the  points.  Central-difference  formulas  are  used  to  approximate  the 
term,  and  this  approximation  leads  to  the  following  system  of 
algebraic  equations,  whose  solutions  approximate  the  solution  of  problem 
(2.2)  at  the  set  of  discrete  points  x^  ; 

A(|)  ■ f . (2.5) 

•v 

Here,  A Is  a trldlagonal  matrix  of  order  J , 


*Least-squares  methods  for  the  solutions  of  Indefinite  equations  may 
accelerate  the  convergence  significantly. 
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A - [1  -S  l]  , 

2 2 

S - 2 - , 

• 

^ is  ths  vsctor 

whose  J^‘*  element  4)^  approximates  the  solution  of  problem  (2.2)  at 
the  point 


! - [-^1 


where  £ denotes  the  transpose  of  £ and  where 


Xj  ■ jAx  . 

The  system  of  equations  that  results  from  multiplying  eq.  (2.5)  by  A*  Is 

= A^$  - Af  = fj  . (2.6) 

When  applying  the  Gauss-Seldel  Iterative  method  (point  relaxation)  to 

st^ 

eq.  (2.6),  we  calculate  the  nfl  Iteration  with  the  relation 


where 


Here,  A 


, J - 1,  2,  ...  J , 

^^i-2^T-2  - 

^®^J+l‘*’j+l  ‘‘‘  ‘^j+2*^J+2  ” ’ 

j - 1,  2,  •••  J . 

Is  the  relaxation  factor, 

+ 2 , 2 < J < J-1 

’ 2 

I + 1 , J - 1 , J 
f 1 , 1 < J < J 


J < I , j > J 
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By  numerical  examples  we  demonstrated  that  relaxation,  when 
applied  to  eq.  (2.6),  converges  for  those  values  of  u)  that  cause  It  to 
fall  when  applied  to  eq.  (2.5).  In  each  case  we  calculated  the  direct 
solutions  of  eqs.  (2.5)  and  (2.6).  As  expected,  they  were  Identical. 

We  also  obtained  the  exact  solution  of  problem  (2.2),  and,  by  comparing 
It  to  the  solution  of  eqs.  (2.5)  or  (2.6),  we  were  able  to  determine  the 
effect  of  the  truncation  error.  We  applied  point  relaxation  with  - 0 , 

(t)j^  * 1 , and  J ■ 19  to  both  systems  (2.5)  and  (2.6).  For  a value  of 
u)  below  IT  (o)  > Tr/2)  , we  found  that  both  cases  converge.  For  a value 
of  (I)  greater  than  ir  (u  - 3ir/2)  , we  fo\ind  that  the  iterative  method 
converges  when  applied  to  system  (2.6),  but  diverges  when  applied  to 
system  (2.5).  Table  I shows  the  results  obtained  for  thla  case  (u  ■ 3Tr/2) . 
The  first  column  of  the  table  gives  the  value  of  x^  at  which  the  solution 
is  found.  The  second  column  shows  the  point-relaxation  (Gauss-Seldel) 
solution  of  eq.  (2.6).  Convergence  was  slow;  It  was  achieved  after  2937 
Iterations.  Convergence  Is  attained  when  the  residual 

<l>?  1 - 2(^"  + ()>“ 

max  |— ^ u 4)"|  < 0.01  , J - 1,  2,  •••,  J . 

j bai  J 

The  relaxation  factor  X Is  equal  to  Its  optimum  value  1.82.  The  third 
column  of  table  I gives  the  solution  of  both  systems  (2.5)  and  (2.6)  obtained 
by  Gauss  elimination,  and  the  last  column  gives  the  exact  solution  of 
the  boundary-value  problem  (2.2).  The  difference  between  the  values  In 
the  last  two  columns  Is  attributable  to  the  truncation  error;  l.e..  It 
Is  attributable  to  approximating  the  differential  equation  by  a set  of 
algebraic  equations.  As  the  value  of  u decreases,  the  truncation 
error  decreases,  and  the  rate  of  convergence  Increases.  For 
0)  - ir/2  the  difference  between  the  exact  solution  and  the  solution  of 
the  approximating  matrix  equation  Is  less  than  0.02Z.  Moreover,  conver- 
gence was  achieved  after  1629  iterations  when  equation  (2.6)  was  solved 
by  relaxation  methods,  and  It  was  achieved  after  44  iterations  when 
equation  (2.5)  was  solved  by  relaxation  methods.  This  result  Indicates  the 
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slowness  of  the  A*A  method  In  comparison  to  the  usual  relaxation 
method  when  It  converges.  In  all  the  examples,  the  optimal  values  of 
the  relaxation  factors  were  used. 

The  above  method  can  be  extended  easily  to  two-dimensional 
problems.  The  slowness  of  Its  convergence,  however,  makes  It  nonpractlcal. 
An  alternative  relaxation  procedure  for  solving  Helmholtz's  equation  for 
frequencies  beyond  the  critical  frequency  Is  given  below. 

2.2.2  The  Bisecting  Method 

For  given  values  of  u)  , Ax  , and  Ay  , there  Is  a family  of 


critical  values  L, 


vcAuco  * ^2cr  * values  of  and  L12  outside 

this  critical  range,  relaxation  procedures  diverge.  Therefore,  If  we 
wish  to  solve  equation  (2.1)  In  a rectangle  having  dimensions  beyond  the 
critical  range,  we  may  divide  the  rectangle  Into  smaller  rectangles,  each 
having  dimensions  below  the  critical  values.  For  given  values  of  ^ on 
the  boundaries  of  the  rectangles.  Iteration  by  relaxation  Is  certain  to 
converge  In  each  of  them. 

To  solve  the  original  problem  (2.1)  we  are  required  to  match 
the  solutions  In  the  different  subregions  at  the  Internal  boundaries  where 
overlap  occurs.  These  matches  are  achieved  by  Iterating  upon  the  boundary 
conditions  of  the  subregions  (the  outer  Iteration).  The  Inner  Iterations 
for  each  subregion  are  done  by  the  usual  relaxation  method. 

The  above  Idea  of  dividing  the  original  region  Into  subregions 
for  values  of  o)  greater  than  can  be  applied  to  the  one-dlmenslonal 

problem  (2.2).  Equation  (2.2a)  Is  approximated  by  a set  of  algebraic 
equations  at  J discrete  points  Xj(-jAx)  , j - 1,  2,  *••  , J In  the 
Interval  of  Interest  P . Th«  original  Interval  Is  divided  Into  two 
subintervals  and  P2  , each  having  a length  below  the  critical 
length,  where 


^ ^ ° ^ ^ *JMAX1  + 1 

^2  ‘ *JM1N2  - 1 1 * 1 ^ » 


JMIM2  - JMAXl  4-  1 . 
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For  the  n outer  iteration,  the  boundary  conditions  Impoaed  on  the 
aolution  for  Interval  are 


♦(0)  - 


(2.7a) 


“ ’^1 


(2.7b) 


The  corresponding  boundary  values  Imposed  on  the  solution  for  interval 


*JMIN2-1^  " ’•'2 

<I>(1)  - <Pr,  . 


(2.8a) 


(2.8b) 


In  Interval  Pj^  M relaxation  sweeps  are  made  towards  solving  equation 
(2.2a)  with  boundary  conditions  (2.7).  We  then  make  M relaxation  sweeps 
for  solving  the  same  equation  in  Interval  P2  • with  boundary  conditions 
(2.8)  Imposed.  The  optimal  relaxation  factors  are  used  in  both  cases. 
The  Iterative  solution  at  point  J at  the  m^^  inner  iteration  and  the 
n*"^  outer  iteration  is  denoted  by  • Th*  values  for  the  inner 

boundary  conditions  and  <|«2  are  given  below.  They  are  chosen  to 

give  the  exact  solution  to  our  problem  in  two  outer  Iteration  steps, 
provided  the  inner  iterations  produce  an  exact  solution  to  the  problems 
in  the  subregions  P^  and  P2  . The  linearity  of  the  problem  is  there- 
fore exploited: 


where  and  C2  are  arbitrary  unequal  constants. 


+ (1  - 


n - 3,  4,  , N 


r 


, N 


and 


a“  - 


.M,n  ,n 
**’JMIN2  ■ '*'l 


(A^»“  _ W;“)  _ 


The  solution  after  the  n^^  outer  Iteration  Is  given  by  4°  , where 

*]~  *]  » J - 1.  2,  •••  , J 

.n  n.M,n-l  . ,,  n» ,M,n 

♦j  - a (j>j’  + (1  - a 

J - 1,  2,  •••  , J 
n • 2,  3,  •••  , N 

In  our  numerical  example  we  have  made  the  following  specific 

choices 


♦l-0 

♦r-  1 

J - 19 
JMAU  - 10 
JMIN2  • 11 
Ax  ■■  0.05 


C^-l 


The  problem  was  solved  for 


and  the  method  was  found  to  converge  in  a reasonable  nuid>er  of  Iterations. 
We  assumed  convergence  when  the  residue  was  less  than  0.01,  where  the 
residue  after  the  n^**  outer  iteration  is  given  by 
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Table  II  shows  the  effect  of  changing  M on  the  rate  of  con- 
vergence. To  compare  the  rate  of  convergence  of  this  method  with  the 
usual  relaxation  method,  we  calculated  an  example  with  u > ir/2  . The 
results  are  given  In  table  III.  Convergence  for  the  usual  relaxation 
method  occurs  In  44  Iterations.  The  present  method  converges  In  only 
40  Iterations,  when  the  number  of  Inner  Iterations  per  outer  Iteration 
Is  5. 

The  bisecting  method  therefore  seems  to  be  suitable  for  solving 
the  problem  of  Interest  at  frequencies  beyond  the  critical  value.  The 
fact  that  the  rate  of  convergence  Increases  by  reducing  the  number  of 
Inner  Iterations  per  outer  Iterations  and  the  fast  convergence  desion- 
strated  by  the  one-dlmenslonal  example  leads  us  to  believe  In  the  possibility 
of  the  extension  of  this  method  to  two-dimensional  problems  with  no 
serious  obstacles. 

2.2.3  The  Transform  Method 

We  mention  briefly  here  that  J.  R.  Molberg  (1968)  developed.  In 
his  thesis,  a transform  technique  which  yields  solutions  to  problems  for 
which  the  Iterative  process  Is  unstable.  He  used  this  method  successfully 
to  solve  problems  (2.1)  and  (2.2). 

2.3  Other  Methods  for  Solving  Helmholtz's  Equation 

Ue  have  described  above  two  modified  relaxation  procedures  for 
solving  the  Helmholtz  equation  for  frequencies  beyond  the  critical 
frequency.  These  methods  were  suggested  as  remedies  for  the  divergence 
problem.  Other  methods  that  can  be  used  to  solve  Helmholtz's  equation 
are  described  below. 

2.3.1  The  Shooting  Method 

Let  us  consider  the  one-dlmenslonal  problem  (2.2).  To  solve 
this  problem  by  the  shooting  method,  we  obtain  the  solutions  and  , 

each  of  which  satisfies  equation  (2.2a)  and  Initial  condition  (2.2b). 

In  addition,  4)^  and  satisfy  the  Initial  conditions 
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\ * 

where  3^  and  3^,  are  two  arbitrary  unequal  constants.  The  desired 
solution  problem  (2.2)  Is  then  given  by 


where 


(j>(x)  - cu^^(x)  + (1  - a)4)^j(x) 
“ “ •>*(!)  - ♦fed)  ' 


The  numerical  scheme  used  to  solve  eq.  (2.2a)  Is  given  by 


- (2  - u)^AxS(|.^  - , j - 1.  2,  — , J , 

where  ((i^  - 

and  - 2Ax  ^(0) 

The  solution  obtained  for  ^ by  shooting  Is  equal  to  that 
obtained  by  direct  Inversion  except  for  round-off  errors.  It  would 
appear  therefore  that  the  shooting  method,  which  requires  no  Iteration 
for  the  one-dljmenslonal  example.  Is  the  most  attractive  of  all  methods 
suggested.  We  point  out  that  although  the  shooting  method  Is  superior 
to  the  A*A  method  and  the  bisecting  method  for  problem  (2.2),  it  loses 
some  of  Its  attractiveness  In  the  case  of  the  two-dimensional  problem 
(2.1). 

The  two-dimensional  generalisation  of  the  shooting  method  for 
solving  two-point  boundary  value  problems  is  only  useful  for  small 
arrays  because,  in  this  method  It  Is  necessary  to  Invert  a badly  ill- 
conditioned  matrix,  and  the  degree  of  111-condltloning  Increaaea  with 
the  array  else  (see  McAvaney  et  el.  (1971)).  A method  which  overcomes 
this  difficulty  was  devised  by  Dietrich  et  al.  (1975).  It  consists  of 
dividing  the  region  of  computation  Into  a nux^er  of  small  subregions. 
Certain  boundary  conditions  are  aasuaed  on  the  interior  boundaries,  and 
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solutions  are  obtained  In  each  subregion  by  using  the  generalized  shooting 
method.  An  Iterative  process  Is  then  defined,  and  the  values  at  the 
Interior  boundaries  are  modified  after  each  Iteration. 

2.3,2  The  Direct  Method 

We  note  that  the  finite-difference  approximation  to  equation 
(2.1)  may  be  written  In  the  form 


IXj)  - f , (2.9) 

where  D Is  the  coefficient  matrix  and  a block  trldlagonal  matrix.  The 
superdlagonal  and  the  subdiagonal  blocks  are  the  Identity  matrix,  and 
the  diagonal  blocks  are  trldlagonal  matrices.  A direct  solution  of 
equation  (2.9)  may  be  obtained  easily  by  the  method  described  by 
Isaacson  and  Keller  (1966). 

2.3.3  Extrapolation  Method 

This  method  avoids  the  unstablllty  problems  by  solving  eq. 

(2.1)  for  values  of  to  below  the  critical  frequency.  Approximate 
solutions  for  values  of  to  greater  than  the  critical  frequency  are 
obtained  by  extrapolation.  A description  of  this  method  Is  given  In 
appendix  A. 

2.4  Concluding  Remarks 

Methods  for  solving  the  unsteady  harmonic  small-perturbatlon 
transonic  flow  equation  have  been  discussed.  Remedies  for  the  divergence 
problem  at  frequencies  beyond  Che  critical  frequency  have  been  proposed, 
and  these  are  shown  to  work  by  testing  them  on  one-dlmenslonal  examples. 
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3.  SHOCK  MOVEMENT  FOR  THE  HARMONIC  APPROACH 

3.1  Introduction 

In  this  part,  we  Investigate  the  problem  of  the  shock  movement 
caused  by  unsteady  perturbations.  In  Weatherin' s and  Traci's  work,  the 
shock  Is  not  perturbed  from  the  position  calculated  from  the  steady 
solution.  For  some  problems,  however,  the  shock  motion  Is  Important  and 
affects  the  whole  flow  field.  To  account  for  these  effects,  we  derive, 
from  a perturbation  expansion  for  the  shock- jump  conditions,  an  equation 
representing  the  shock  movement.  The  appropriate  jumps  at  the  unper- 
turbed (steady)  shock  position  are  obtained  by  an  analytical  continuation 
of  upstream  and  downstream  conditions.  To  account  for  the  shock  movement, 
we  Impose  the  appropriate  jumps  on  the  perturbation  potential  with  a shock- 
fitting procedure.  A two-dimensional  example  Is  computed  to  Illustrate 
the  procedure. 

For  steady  perturbations  calculated  by  an  Integral  equation 
method,  Nixon  (1977)  has  reported  an  alternative  approach  using  the 
method  of  strained  coordinates.  In  this  paper  we  adopt  the  same  ap- 
proach for  small  harmonic  perturbations  computed  by  finite-difference 
methods.  In  this  approach,  the  coordinates  are  strained  In  such  a way 
that  the  shock  Is  always  fixed  at  Its  steady-state  position.  The  per- 
turbation of  the  jump  conditions  (In  the  strained  coordinates)  yields 
the  same  equations  as  the  weak  solution  of  the  (linear)  perturbed  equa- 
tion; hence,  no  shock  fitting  Is  necessary.  The  perturbed  equation, 
however.  Is  more  complicated  because  It  has  nonhomogeneous  terms  ac- 
counting for  the  shock  movement.  Also,  the  boundary  condition  for  the 
perturbed  potential  Is  altered  (the  airfoil  Is  distorted  so  that  the 
shock  location  Is  unchanged  by  the  perturbation). 

Finally,  we  show  that  the  strained-coordinate  method  Is 
equivalent  to  the  direct  method  of  transferring  the  jump  conditions  to 
the  steady-state  shock  location. 

3.2  Basic  Governing  Equations 

The  unsteady  transonic  small-dlsturbance  equation  can  be 
written  In  the  form 

<»t  + - W.  - T ^ <Vy  • ”•*> 
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where  the  term  Is  neglected  for  low  frequencies.  The  jump  con- 

ditions, admitted  by  the  weak  solution  of  eq.  (3.1),  are 

(3.2a) 
(3.2b) 


where  <(^>  and  |(|)|  denote  the  average  and  the  jump  of  <p  across  the 
shock  X - X°(y;t).  To  complete  the  formulation  of  the  problem,  we  must 
Include  the  tangency  boundary  condition,  a far-fleld  behavior,  and  a 
Kutta  condition  for  lifting  airfoils. 


For  small  harmonic  perturbations,  we  let* 


(j»  ■ (J)  (x,y)  + eRe 


e (p  (x,y) 


+ 


2 „ 
c Re 


/■a2.0  . 21a)t  .2,2. 

(♦  ’ + e ) 


(3.3) 


where  the  oscillating  airfoil  has  the  following  boundary  condition: 


<l>y(x,0)  - ^^(x)  + eRe  |e^^fj^(x)  j . (3.4) 

The  zeroth-order  problem  Is  given  by 

(K  - - 0 , (3.5.) 

♦°<x,0)  - , (3.5b) 

and  the  first-order  problem  (where  the  6 term  Is  neglected)  Is 

(K  - ♦°)*^  + (3.6.) 

(j.J(x,0)  - fj(x)  . (3.6b) 


*Thls  form  Is  equivalent  to  ^ ■ 


+ ♦^.O*  ^ ^2,2^21a,t  ^ ^2,2*^-21a)t 


). 
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The  boundary  conditions  (3.Sb)  and  (3.6b)  are  applied  In  the  airfoil  mean 
surface.  Note  that  eq.  (3.6a)  Is  linear  with  discontinuous  coefficients. 
Its  solution  admits  jumps  only  at  that  place  where  the  steady  shocks 
occur.  We  discuss  this  problem  below. 


Equation  for  Shock  Movement 


3.3.1  Derivation  of  the  Equation 

The  problem  can  be  examined  In  two  ways.  First,  we  perturb 
the  shock  and  transfer  the  perturbed  jump  conditions  to  the  old  position 
(the  steady  state).  These  transferred  conditions  are  Imposed  on  the 
perturbation  potential. 

In  the  second  approach,  following  Nixon  (1977),  the  method  of 
strained  coordinates  Is  used.  The  coordinates  are  strained  so  that  the 
shock  Is  always  fixed  at  Its  steady-state  location  (In  the  strained 
coordinates).  The  differential  equation  governing  the  perturbation 
potential.  In  the  strained  coordinates,  contains  nonhomogeneous  terms 
that  account  for  the  shock  movement.  Also,  the  boundary  condition  Is 
distorted.  This  approach  Is  presented  In  section  3.4. 

Following  Cheng  and  Hafez  (1973)'*,  the  shock  Is  perturbed  In 
the  same  way  the  flow  field  Is  perturbed;  namely,  let 


X°(y;t)  - X°(y)  + eRe[e^V(y)] 


+ e Re  (X 


2,0  ^ ^21^tx2,2  ^ 


(3.7) 


The  jump  conditions  (eqs.  (3.2a)  and  (3.2b))  are  given  In  terms  of  the 
averages  and  the  jumps  across  the  unknown  perturbed  shock  X^(y;t)  . 
The  appropriate  jumps  at  the  unperturbed  shock  X°(y)  may  be  obtained 
by  analytical  continuation  of  upstream  and  downstream  conditions.  For 
example,  consider  the  jump  In  a quantity  U at  X^  , namely  |]u]|  ^ . 
Up  to  first  order,  this  jump  can  be  expressed  In  terms  of  |u]|  ^o  , 
l^xOxO  shown  In  figure  1. 


*M.  Hafez  would  like  to  thank  Prof.  H.  K.  Cheng  of  U.S.C.  for  an  Interesting 
discussion.  A similar  approach  Is  studied  In  Cheng  and  Hafez  (1973),  where 
It  Is  applied  to  a lift  perturbation  problem. 
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Figure  1.  Perturbation  of  a Diacontinuoua  Function 


We  can  write 


JiiljjD  - Oj  + (ii„  + + ... ) 

■ «d- “«>  + «*  K)d-<v»h  — 

- IdIjO  + ex*  0»J  „ + •..  . 

X 

similarly,  we  have  <U>  _ ■ <U>  + eX^<U  > + ••• 

x°  x°  * x° 

Now,  consider  condition  (3. 2b) expanded  in  terms  of  e : 

|[(j)°  + ecj)^  + *’*D^  “ ’**1  o 


+ eX^lI^°  + e<l>i  + ‘“I  yo  + •••  - 0 . (3.8) 

XX  A 

The  zeroth  and  first-order  relations  are  given  by, 

for  e°,  - 0 . (3-9) 

for  l[(^^DjjO  - -X^  (l^°JjjO  . (3.10) 

Similarly,  condition  (3.2a)  gives, 

j 

tore”,  <K  - - -(1^)'  . (3.11) 

for  e‘,  210.X'  - + x‘4  - 2 ^ (3.12) 


Equation  (3.12)  is  an  ordinary  differential  equation  in  X^  . Together 
%d.th  eq.  (3.10),  it  imposes  the  proper  Jump  conditions  on  . A shock- 
fitting procedure  is  necessary  since  the  weak  solution  of  the  (linear) 
perturbed  system  (3.6)  admits  discontinuities  other  than  the  conditions 
that  result  from  perturbing  the  Jumps  admitted  by  the  fully  nonlinear 
system.  (3.1)  To  calculate  the  pressures  in  the  region  of  the  shock,  we 
must  take  into  account  the  motion  of  the  shock  as  discussed  in  appendix  B. 


3.3.2  Slapllflcatlon  for  Normal  Shock*  ^ 

For  a locally  normal  shock,  the  term  ®‘l*  (3.12)  is 


neglected;  hence, 

21ajX^  - ^4)^  + 4)°  \ 

\^x  xx/^o 

Eliminating  between  (3.10)  and  (3.13),  we  have 

KB  .0  *<21“  - *lxlo  0*'1I  .0  ■ “ • 


(3.13) 


(3.14) 


Equation  (3.14)  is  the  jump  condition  to  be  imposed  on  . After 
is  calculated,  we  can  evaluate  X^  from  sq.  (3.13)  or  (3.10). 


3.3.3  A Shock-Fitting  Procedure  for  the  4>^  Problem 

First,  we  solve  the  steady  problem  , using,  for  example, 
Murman's  fully  conservative  schemes.  From  the  4>°  solution,  we  let 
|[4)°J  jjO  “ a and  ^21i»)  - ^ » where  X°  is  considered  part  of 

the  4>°  solution.  The  first  point  downstream  of  the  shock  X°  is 
identified  by  S.P.  in  figure  2.  At  this  point,  we  replace  the  finite- 
difference  equation  for  4>^  by 
the  following  relation: 


+ b(4)J  - 4»i.i)  - 0 . 


(3.15) 


Equation  (3.15)  is  a first-order  approximation  of  eq.  (3.14). 

The  shock  is  treated  as  an  internal  boundary  with  a boundary 
condition  in  the  form  of 

a/  + bj**  - Cj  . 

where  a^  , b^  , and  c^  are  constants,  (+)  denotes  the  conditions  just 
downstream  of  the  shock,  and  c^  contains  the  conditions  on  the  upstream 
side  of  the  shock. 


*For  normal  shock,  ^4)^^-  0 . 
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Figure  2.  Shock  Point*  and  Shock  Location 
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3.3.4  Computatlohal  Example 

Using  Che  above  theory,  we  can  compute  a two-dimensional  example 


of  a pulsating  parabolic  arc  airfoil.  The  airfoil  Is  placed  on  Che  line 
y - 0 , and  Its  upper  surface  Is  governed  by  the  equation 

y^  ■ 0.2  X (1  - x)(l  +0.1  sin  a>t)  , 

where  o)  ■ 0.818  . 

We  now  solve  the  problem  In  the  rectangular  domain 

-1  < X < 2 , 

0 < y < 2 . 

We  set  the  freestream  Mach  number  equal  to  0.818  and  solve  the  lowest 

order  steady-state  problem 

[(1  - m2)  - (Y  + ° 

♦°(-l.y)  - 0 . ♦°(l,y)  - 0 

♦°(x,2)  - 0 

0 X < 0 

(<>y(x,0)  - 0 X > 1 

lo.2(l-2x)  0 < X < 1 


by  using  the  Murman  type-dependent  relaxation  scheme.  A plot  of  Che 
solution  Is  shown  In  figure  3.  The  unsteady  perturbation  equation 


Fifurc  3.  Steady-State  Solution  for  a Parabolic  Airfoil 


with  the  boundary  conditlona 


♦^(-l.y)  - 0 ♦^l.y)  - 0 

♦^x,2)  - 0 


0 X < 0 

♦y(x,0)  - 0 X > 1 

.02x(l  - x)  0 < X < 1 


can  then  be  solved  by  using  central  differencing  in  subsonic  regions  and 
backward  differencing  in  supersonic  regions.  At  the  shock  point,  the 
double  derivative  appearing  in  eq.  (3.16)  can  be  expressed  as 


- i- 

^xx  6x 

We  then  approximate  the  term 
it  equal  to 


^i+1  - ^i  .+ 

S ♦x  • 

2 2 
, »Ll  - < 

*^x  ^ Ax 


from  eq. 


(3.15)  and  set 


^i-1 


- ♦ 


i-2 


Ax 


♦i-i> 


This  is  necessary  to  satisfy  both  the  boundary  condition  and  the  shock 
jump  relation  at  the  base  of  the  shock. 

Figure  4 shows  the  shock  position  obtained  from  the  present 
calculation.  Also  shown  in  the  figure  is  the  shock  position  obtained  by 
solving  the  unsteady  small-disturbance  equation  (Ballhaus) . The  shift 
in  the  shock  mean  position  is  caused  by  higher  order  terms  which  are 
unaccounted  for  in  this  calculation.* 

3.4  An  Alternative  Approach  Using  the  Method  of  Strained  Coordinates 

3.4.1  Formulation 

Instead  of  eq.  (3.3),  we  use  an  expansion  of  in  strained 
coordinates  (s,T)  as  follows: 


*A  modified  quasi-steady  approach  may  be  used  where  the  unsteady  (^  ) term 

is  evaluated  by  using  the  first-order  perturbation  solution.  This  Spproach 
leads  to  a nonhomogeneous  equation  of  the  form 

(K  - ((»  )())  +4  - (Re  2ia)M^())^e^“^ 

^x  ^xx  ^yy  “^x 
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Yy  = 0.2Xn  - x)(1  + 0.1  sin  a>t) 


Pulsating  Airfoil 


<l>(x,y;T)  - <|)^(8,y)  + eRe  | e^4>j^(8,y)  j + •••  , 
X - 8 + eRe  je^Xj^(8)j  + •••  , 


t - T . 


(3.17) 


(3.18) 


(3.19) 


In  particular,  the  ahock  location  in  the  physical  plane  is  related 

to  the  shock  location  in  the  strained  coordinates  by  the  following 


relation: 


X°  - S°  + eRe  e^‘^Xj^(S°)]  + •••  . 


(3.20) 


That  is,  X^(S  ) is  the  shock  movement. 


The  transformation  derivatives  are  given  by 


1^-  (1  - eRe  * 

It  " I?  + •••  >|i-  • 


(3.21) 


(3.22) 


Using  (3.17)  and  (3.18)  in  eq.  (3.1)  (with  K 1,  and  neglecting  the  8 
term)  yields 


(l-(j)  )(^  +0  “0, 

os  OSS  oyy 

and  (j),  4’i  “ (4’  ) + 2i(i)(4ii  ~ X,  4 ) 

^Iss  ^lyy  ^os^ls  s ^1  l^os  s 

+ [x-  (4  - 4^  )1  + X,  (4  - >s4^  ) , 

I la'^os  ^os'Js  Is'^os  ^^os'^s  ’ 

while  using  (3.17)  and  (3.18)  in  eqs.  (3.2a)  and  (3.2b)  yields 


<4„.>  - 1 , 
08 


J-o  . 


<'>l8>  - Xi,(S°)  + 21(UX^(S°)  , l4iE  - 0 . 


(3.23) 


(3.24) 


(3.25) 


(3.26) 
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Note  that  eqs.  (3.25)  and  (3.26)  are  the  jump  conditions  admitted  by  the 
weak  solutions  of  eqs.  (3.23)  and  (3.24),  respectively.  The  boundary  conditions 
In  the  physical  plane  are  transformed  to 

- W^(o)  on  the  airfoil,  and  (3.27) 

♦ly(s,0)  ■ Wj^(e)  + Xj^(8)W^^(o)  on  the  airfoil  . (3.28) 

Equations  (3.23),  (3.25),  and  (3.27)  determine  the  steady  state,  while  eqs. 

(3.24),  (3.26),  and  (3.28)  determine  the  perturbation  potential  In  the  strained 
coordinates . 

3.4.2  Difference  Equations 

For  the  steady-state  problem,  centered  differences  are  used  In 

the  subsonic  region,  backward  differences  In  the  supersonic  region,  a 

parabollc-polnt  operator  at  the  sonic  line  (to  exclude  expansion  shock) , 

and  a shock-point  operator  at  the  shock  point  to  Impose  the  right  Jump 

condition  (conservation  of  mass) . 

For  the  perturbation  potential,  a parabollc-polnt  operator  and 

a shock-point  operator  are  needed.  Also,  X,  (S  ) Is  unknown. 

1 o 

We  first  consider  Che  parabollc-polnt  operator.  In  (x,y,c) 
coordinates,  ~ ^ X the  sonic  line.  To  guarantee  finite  acceleration 
at  the  sonic  line,  eq.  (3.1)  reduces  to 

V - • (3.29) 

(3.30) 

(3.31) 

In  eqs.  (3.30)  and  (3.31)  centered  differences  are  used  everywhere. 

Equation  (3.24)  at  (♦^^  - 1)  reads 

“ - • ♦...  [♦!.  - «1.  + ) 

! 

... 


Parabollc-polnt  operators  In  (s,y,T)  coordinates  are  given  by 


V ■ • 


To  keep  finite  at  the  (4>  - 1),  we  must  satisfy  the  following 

X98  08 

condition: 


♦is  “ *1.  + . at  <0^^  - 1 . 


(3.33) 


Equation  (3.33)  Is  used  to  determine  the  shock  movement  A (assuming,  for 
example. 


Xj^(s)  - Asd  - 8)/S°(l  - S°) 

as  In  Nixon's  work  (1977)).  Note  that  eq.  (3.33)  la  consistent  with  eq. 
(3.26)  when  the  shock  strength  vanishes. 

We  next  consider  the  shock-point  operator.  Equation  (3.24)  Is 
written  In  a conservative  form.  A fully  conservative  scheme  that  admits 
the  appropriate  jumps  eq.  (3.26)  Is  written  as  follows: 


(Shock-Point  Operation)^  ■ (Elliptic  Operator) + (Hyperbolic  Operator)^  . 

(3.34) 


The  y terms  sre  central  differenced. 

After  algebraic  manipulation,  eq.  (3.34)  reduces  to  (neglecting 
y terms  for  a locally  normal  shock) 


where 


^ 2iadCj^(S°)  , 

+ _ ♦1,1+1  ~ ^1,1 

^ AX 

u-  - ♦1.1-1  ~ ♦1.1-2 

AV  . 


0.35) 


Equation  (3.35)  Is  a first-order  approximation  of  Jump  condition 

(3.26). 


3.5  Equivalence  of  the  Two  Methods 

Finally,  we  want  to  show  that  these  results  are  consistent 
with  the  previous  approach  and  that  the  two  approaches  are  equivalent. 
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We  substitute  eq.  (3.18)  into  eq.  (3.3),  together  with  a Taylor  series 
expansion  for  around  the  point  x,y  , and  collect  terms  of  equal 
orders.  The  result  Is 


♦ (x.y.t)  - + cRe  U^s.y)  + Xj(s)$°  (a.y)]  | . 


Comparing  eqs.  (3.36)  and  (3.17),  we  find 


(3.36) 


♦“(s.y)  ■ 4'Q(s,y)  , 

♦^(8,y)  + Xj^(s)4»°(s,y)  - (Ji^(s,y)  . (3.37) 

The  equation  for  ♦^(s,y)  reads 


(3.38) 


with  the  boundary  condition 


♦y(s,0)  - W^(s)  (3.39) 

on  the  airfoil  and  the  shock-jump  conditions 

- -Xj(S°)  , (3.40) 

+ Xi(S°)<;4.°g)  - 2ia>X^(S°)  . (3.41) 

Equations  (3.38),  (3.39),  (3.40),  and  (3.41)  were  used  before  with  a shock- 
fitting procedure  to  determine  . Hence,  the  two  systems  are  equivalent. 

3.6  Conclusions 

We  have  derived  an  expression  for  the  perturbed  shock  position 
for  unsteady  transonic  calculations.  This  result  provides  us  with  a 
consistent  method  for  predicting  the  shock  motion  of  small  unsteady 
harmonic  perturbations.  A two-dimensional  example  has  been  calculated. 
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FINITE  ELEHENT  METHODS  FOR  STEADY  AND  UNSTEADY  TRANSONIC  FLOWS 


4.1  Introduction 

For  Isentroplc,  Irrotatlonal,  transonic  flows,  Euler's  equations 
reduce  to  the  mass  conservation  equation.  Evaluating  the  density  from 
the  Bernoulli  (energy)  equation,  we  find  the  flow  Is  governed  by  a single 
second-order  nonlinear  hyperbolic  equation  In  terms  of  a velocity  potential. 
The  steady  equation,  however,  Is  of  mixed  type  (l.e.,  locally  hyperbolic 
(elliptic)  for  the  supersonic  (subsonic)  region.  The  change  of  the  type 
of  the  equation  and  the  existence  of  a dlscontlnous  solution  (shock 
waves)  are  attributable  to  the  nonlinearity) . It  seems  that  the  unsteady 
problem  Is  more  amenable  to  analysis.  In  fact,  the  steady  solution  Is 
obtained  as  the  steady-state  limit  (for  large  times)  of  the  unsteady 
problem. 

The  controlling  parameters,  In  general,  are  the  Mach  number, 
geometry  parameters  (angle  of  attack,  thickness  ratio,  aspect  ratio, 
etc.)  and  the  unsteady  disturbance  parameters  (frequency  and  amplitude). 

With  standard  limit  processes,  different  approximations  can  be  obtained. 

The  most  famous  of  these  Is  the  transonic  small-disturbance  theory.  For 
the  high-frequency  range,  the  nonlinear  terms  can  be  neglected,  and  the 
equation  Is  linear.  For  the  low-frequency  range,  the  nonlinear  term  is 
Important,  and  other  linear  terms  like  tUa  term  are  not.  Although 

the  approximations  lead  to  great  simplification,  some  mathematical 
difficulties  are  associated  with  the  resulting  governing  equation  In 
this  range  (a  parahyperbollc  equation  with  one  real  characteristic) . 

The  solution  procedures  for  the  unsteady  problems  are  numerous. 

On  one  hand,  the  straightforward,  quasi-steady  approach,  where  time  is 
considered  a parameter,  may  be  Justifiable  In  the  low-frequency  range. 

If  we  have  an  efficient  method  for  solving  the  steady  problem.  On  the 
other  hand,  direct  time  Integration  of  the  unsteady  equation  allows  for 
the  full  unsteady  effects  to  take  place.  In  between  these  two  methods 
are  the  linear  perturbation  methods  for  small  unsteady  disturbances, 
namely,  the  frequency  and  Indlclal  response  methods.  In  the  frequency 
decomposition  method,  a sequence  of  steady-llke  problems  (for  different 
frequencies)  are  solved;  while,  for  the  Indlclal  method,  the  response  of 
a linearized  time-dependent  equation  to  an  Impulse  Is  calculated.  As  men- 
tioned earlier  In  this  report,  perturbation  of  a system  with  a discontinuous 
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solution  must  be  treated  with  a special  consideration.  Nonlinear  finite- 
amplitude  effects  can  be  Included  If  the  main  flow  Is  allowed  to  Interact 
with  the  disturbance  leading  to  a coupled  system  of  nonlinear  equations 
for  the  main  flow  as  well  as  the  perturbations. 

Numerical  methods  are  needed  to  Implement  these  procedures. 

In  the  following  sections,  application  of  the  finite-element  method  to 
the  steady  and  unsteady  transonic  problems  is  studied. 

In  recent  years  various  finite-element  techniques  have  been 
developed  for  steady  small-disturbance  transonic  flow  problems.  Chan, 
Brashears,  and  Young  (1975)  presented  the  first  application  of  finite- 
element  methods  In  transonic  flow.  They  resorted  to  a least-squares 
model  after  discovering  that  Galerkln  formulations.  In  conjunction  with 
an  Iterative  scheme  based  on  a sequence  of  Poisson  solutions,  did  not 
converge  for  supercritical  flows.  Wellford  and  Hafez  (1976)  constructed 
a workable  Galerkln  model  by  introducing  a mixed  variational  principle 
and  an  appropriate  direct-gradient  solution  algorithm.  Hafez,  Murman, 
and  Wellford  (1976)  developed  a Galerkln  model  that  accounted  for  the 
change  In  type  of  the  differential  equation  In  a manner  similar  to  the 
technique  used  in  Murman* s original  work.  Glowlnskl,  Perlaux,  and 
Flronneau  (1976)  Introduced  an  optimal  control  formulation  for  transonic 
flow.  The  finite-element  method  has  not  been  used  for  unsteady  transonic 
flows  however,  there  is  some  work  in  the  literature  on  the  solution  of 
shallow  water  equations,  which  are  similar  to  the  transonic  equation 
(Y»2). 

In  the  present  study,  we  have  constructed  an  alternate  Galerkln 
model  for  the  small-disturbance  transonic  problem  and  have  developed  a 
solution  algorithm  for  the  resulting  equations.  The  problem  Is  formu- 
lated In  terms  of  velocities,  and  an  Implicit  finite  element  Is  constructed 
by  using  the  time-dependent  model  of  Magnus  and  Yoshlhara  (1970)  for  both 
steady  and  unsteady  problems.  Because  the  resulting  algorithm  Is  essen- 
tially independent  of  the  type  of  the  flow  (subsonic  or  supersonic),  it 
is  well  suited  for  finite-element  models.  Curved  boundary  shapes  and 
Irregular  meshes  can  be  used.  The  advantages  of  the  finite-element 
method  In  modeling  highly  Irregular  domains  and  in  incorporating  higher 
order  accurate  elements  can  therefore  be  utilized. 
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In  the  following  sections  the  physical  problem  Is  posed,  the 
Iterative  solution  algorithms  are  Introduced,  the  convergence  criteria 
for  the  solution  algorithm  are  determined,  and  numerical  results  are 
presented . 


4.2  The  Physical  Problems 

We  consider  the  flow  of  an  Invlscld,  compressible  fluid  about 
a thin  airfoil.  We  adopt  here  the  small-dlsturbance  theory  (see  Landahl 
(1961)).  Let  0 be  an  Infinite  domain  composed  of  points  (x,y)  . Let 
IT  be  the  upper  half  plane  defined  by  y ^ 0 , for  all  x . We  consider 
non-llftlng  symmetric  flows,  and  thus  we  model  domains  S)  such  that 
0 C IT  . 0 Is  the  same  as  (2  except  for  the  region  occupied  by  the 
airfoil.  We  describe  the  boundary  9(2^  of  the  upper  half  plane  of  the 
symmetric  airfoil  by  y > g(x,t)  where  y - 0 Is  the  chord  of  the 
airfoil.  The  boundary  9(22  of  the  flow  Is  the  boundary  (y  ~ 0)  of 
the  half plane  IT  for  points  outside  the  airfoil.  The  boundary  9(2^  of 
the  flow  Is  the  boundary  of  the  half plane  at  <»  In  the  y direction. 

The  boundary  9(2^  of  the  flow  Is  the  boundary  of  the  halfplane  at 
X ■ <*  and  X “ -<*>  . The  boundary  9(2  of  (2  Is  9(2  ■ 9(2j^U  9122^^  9(22U  9(2^ 
Let  4)  be  the  perturbation  velocity  potential,  be  the 

free-stream  Mach  number,  and  y be  the  ratio  of  specific  heats.  Then 
the  problem  Is  described  by  the  following  set  of  equations: 


m2(i+y)4>^  I 

2 'x  j 


X 


2M^(j)  ^ + 
“^xt 


M^(Ji 

®^tt 


In 


(2 


+ (1  + " *^*y  “ 0 O® 


♦ fy  “0  on  9(22 

■ 0 on  9(22L^(2^ 


(4.1) 


For  steady  flows,  suppose  that  the  x and  y velocity  components  are 
defined  by  u ■ (>,x  and  v - ^,y  . Then  (4.1)  can  be  written  In  the 
following  form  (a  system  of  first-order  equations): 
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2 2 


In  f2 


(1  + “^8»x  ■ V ■ 0 


on  3(2, 


V - 0 


3J2^j 


u ■ 0 


on  3(2. 


where  - 1 - IQ 
by  requiring  that 


and  1^2  “ Y)*Q  • Irrotatlonallty  is  imposed 


“’y  " "^’x  “ ° in  (2  . (4.3] 

The  low-frequency  unsteady  equation  (omitting  the  4)^^  term)  is  a 
system  of  first-order  equations  in  the  form  of: 

*^2  2 

(Kj^u  - u )x  + Vy  - au^ 

(4.4) 

Uy  - Vx  - 0 

The  high-frequency  equation  (4.1),  can  be  written,  however,  in  the  form 
of  a system  of  first-order  equations  if  we  use  three  variables  u,  v and  w 
where  w ■ ; namely, 

(K^u  - ^ u2)^  + Vy  - 2aw^  - w^ 

V ■■  u 

X t (4.5) 

w ■ v^ 

y t 

2 2 

where  1 - IQ  and  K2  ■ (1+y)M|j^  . Irrotatlonallty  is  imposed 

by  requiring  that 


u - V « 0 
y X 


in  (2 


4.3  Tentative  Solution  Algorithms 

To  construct  a workable  solution  algorithm,  we  introduce  time- 
dependent  and  explicit  artificial  viscosity  terms  as  follows: 


®“*t  " 1 “ 


- *^2  2 

K,u  - -r—  “ j. 
P 2 >x  + 
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+ YU, 


- XU 


14.6) 


XX 


6v,^  ■ u,  - V,  + Xv,  - Fv  . 

»t  »y  »x  ’yy 

0£  course,  the  second-order  terms  In  (4.6)  can  be  associated  with 
viscous  dissipation.  The  terms  Involving  u and  v are  more  difficult 
to  explain  physically.  Lee  and  Kim  (1976)  have  Introduced  such  terms  as 
a correction  to  the  Invlscld  KdV  equation  occurring  In  free-surface 
hydrodynamics.  In  their  theory  the  term  accounts  for  the  no-sllp  boundary 
condition  of  the  viscous  theory.  However,  It  Is  not  necessary  to  have  a 
physical  argument  to  justify  the  perturbation  of  the  original  problem 
Introduced  here.  It  suffices  to  say  that  the  attached  terms  are  mathe- 
matical In  nature  and  produce  a convergent  algorithm  for  supercritical 
flows.  Presumably,  as  the  artificial  viscosity  parameters  are  reduced 
to  zero,  the  correct  approximate  solution  Is  retrieved. 

Introducing  a finite-difference  model  for  the  time  variation  In 
(4.6)  In  association  with  a time  step  At  , we  obtain  the  first  Iterative 
scheme  to  be  studied  here,  which  takes  the  following  form: 


g 

At 


(u 


n+1 


U\ 
- u ) 


„ n4*s  2 n n+lj 
KjU  •*  - y-  u u 


, + v,^*^  + n+^s  n+*i 

* y Yu,^^  - XU  , 


l_(v 

At^ 


n+1 


- V ) - 


n+1 


Xv, 


n-t^s 

yy 


- Fv 


n4^ 


(4.7) 


where,  for  example, 

w * ~ * 

Algorithm  (4.7)  has  a linearization  that  Is  stable  and  second-order 
accurate.  However,  the  nonlinear  algorithm  Is  not  unconditionally 
stable.  Certain  minimum  amounts  of  artificial  viscosity  as  well  as 
restrictions  on  the  time  step  At  are  required  to  stabilize  the  algorithm. 
Note  that  the  use  of  (4.7)  requires  the  solution  of  an  asynmetrlc  linear 
system  at  each  Iteration. 

He  can  obtain  a second  algorithm  by  treating  the  nonlinear 
term  In  (4.6)  differently,  as  follows: 
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X 


(4.8) 


a , n+1  n.  „ n+H  ^2,  n4Ji.2  , n+*s  n4Jj  n+lj 

^(u  - u ) - I^K^u  ^ - ^(u  + v,^  + Yu,j^  - XU  . 


B,  n+1  n.  n+^i  n+*s  , . n+Jj  » n44j 

f-(v  - V ) - u.  / - v,"7  + Xv,^.  - Pv  ^ 


At 


yy 


Algorithm  (4.8)  Is  unconditionally  stable;  however,  it  requires  the 
solution  of  a nonlinear  system  at  each  time  step;  thus,  computationally 
is  not  as  attractive  as  (4.7). 

The  third  algorithm  to  be  considered  is  obtained  by  constructing 
an  explicit  form  of  (4.6).  It  takes  the  following  form: 


a . n+1  n. 
^u  - u ) - 


„ n *^2,  n.2 

Kj^u  - -^-(0  ) 


, + V?  + Yu?  - xu” 

X ’y  ' *xx 


(4.9) 


fl  . n+1  n 


tA\  A A , « A ««  A 

V ) • U.  - V,  + AV,  - Pv 

■y  »X  'yy 


This  algorithm  is  only  conditionally  stable.  In  addition,  certain 
minimum  amounts  of  artificial  viscosity  are  required  to  stabilise  the 
nonlinear  term.  However,  in  some  cases  the  computational  simplicity  of 
this  algorithm  may  be  advantageous. 

The  characteristic  feature  of  these  algorithms  is  the  introduction 
of  explicit  artificial  damping  terms  in  both  space  and  time  for  stability 
and  convengence  considerations.  The  accuracy  of  the  solution  depends, 
of  course,  on  the  slse  of  these  additional  terms.  For  example,  the 
inviscid  solution  can  be  retained  only  if  the  viscous  term  is  of  the 
same  order  as  the  truncation  error.  Similarly,  the  steady-state  aolution 
is  accurate  only  if  the  terms  proportional  to  u and  v , which  %rere 
Introduced  to  dampen  the  Iterative  solution,  are  small.  (These  terms 
can  be  varied  with  Iteration  and  gradually  vanish  when  the  solution 
reaches  the  steady  state.)  The  explicit  introduction  of  these  terms 
facilitates  the  analysis  and  implementation  of  finite  elements.*  Moreover 
the  coefficient  of  the  damping  terms  may  be  chosen  to  optimise  the  rate  of 
convergence . 


For  example.  Instead  of  searching  for  a finite-clement  analog  for  the 
backward  differencing  of  u^  , we  explicitly  add  u to  It  and  use  the 
usual  Galerkln  method  (equivalent  to  centered  differences  for  rectangular 
elements) . 
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For  unsteady  problems,  however,  the  transient  solution  is  of  primary 

concern,  the  same  algorithm  can  be  used  with  proper  choice  of  different 

parameters.  For  example.  If  we  wish  to  simulate  the  transonic  low-frequency 

equation,  0 must  vanish  or  at  least  be  of  the  same  order  as  the  truncation 
2 

error,  namely  (At)  . The  damping  terms  in  time  will  harm  the  solution  and 
therefore  must  be  minimized.  The  main  objective  In  the  choice  of  these 
different  parameters  Is,  In  this  case,  to  minimise  the  phase  error. 

The  explicit  introduction  of  these  "artificial"  terms  enables 
the  algorithm  to  cope  with  completely  different  requirements.  Obviously, 
this  flexibility  Is  advantageous  In  the  study  of  unsteady  transonic  flows. 

4. A The  Weak  Formulation 

Let  L2(n)  be  the  space  of  square  Integrable  functions  on 
n . An  Inner  product  Involving  the  functions  of  L2(0)  Is  defined  as 
follows: 

u,v  ■ uv  dx  dy  . 

2 

The  natural  norm  associated  with  Is  l|u||^  ■ u,u 

Let  !„((})  be  the  space  of  functions  with  the  followi§i  norm: 


In  addition,  let  H (Q)  be  the  Sobolev  space  with  norm 


l•“’xllL2(n)  * ' '“’y' lL2((i)  * 

A weak  form  of  the  system  of  partial  differential  equations 
(4.7)  can  be  defined  as  follows:  Find  u"*^,v"^^  e H^(n)  such  that 


o n+1  n u 

aF  “ ■ “ 


V 2 n n+Jj  „ 

KjU  * - yu  u 


* -X  .'^.w  , 

v"*‘  - y“,S  . (u«*.S)  - 


(4.10) 


+ X v?;^,S  - F v“^,S  , 

for  arbitrary  U,S  e H^(0)  , where  H^(0)  1*  the  space  of  functions  in 
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H (Q)  that  are  zero  on  the  boundary  3R  . Integrating  by  parte  in 
assorted  terms  on  the  right  In  (4.10)  and  using  the  fact  that  W and 
S . zero  on  the  boundary,  we  obtain  the  following  expression: 


- u“,W>  - - 

(4.11) 


- y<v”^.s> 

for  all  W,S  £ H^(fi)  . Equation  (4.11)  represents  the  weak  form  of  (4.7) 
In  a conservative  form,  which  Is  appropriate  for  the  problem  (4.1)  and 
which  naturally  Involves  shocks. 

A weak  form  of  the  system  (4.8)  can  be  developed  using  similar 
techniques.  The  following  form  Is  obtained: 

for  all  W,S  e H^(n)  . A weak  expression  of  the  same  form  can  also  be 
developed  for  (4.9). 

A stability  and  convergence  proof  for  the  Iterative  techniques 
has  been  studied  In  detail  by  Wellford  and  Hafez  (1977). 

4.5  Finite- Element  Models 

In  the  numerical  results  to  be  presented  here,  we  use  algorithm 
(4.11).  To  Implement  this  weak  formulation,  we  subdivide  the  domain 
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E 

n Inco  finite  elements  . Then  Is  U 0 , i<^ere  E Is  the  total 

® ^ 1 ® 
e«l 

number  of  elements  In  the  domain.  On  each  domain  the  velocity 

components  U and  V are  approximated  by  a finite-element  model  of  the 
following  form: 


Ug  - l|)^(x,y)U^  , 

Vg  - lJ»^(x,y)V^  , 


(4.13) 


where  4)^(x,y)  are  finite-element  Interpolation  (shape)  functions,  and 

and  are  the  values  of  u and  v at  the  nodes  of  the  element. 

Of  course,  u and  v could  be  Interpolated  with  different  shape  functions; 

however,  this  procedure  does  not  seem  to  present  any  advantage  In  this 

application.  In  the  numerical  results  to  be  presented  here,  the  eight- 

node  serendipity  element  was  used.  In  particular,  the  mesh  used  In  this 

study  Is  shown  In  figure  5.  While,  In  this  mesh  each  element  has  two 

sides  parallel  to  the  y coordinate  axis,  the  method  demonstrated  here 

should  also  be  valid  for  more  general  mesh  configurations  because  the 

governing  equations  are  hyperbolic  at  all  points  In  the  domain. 

An  approximate  version  of  the  weak  formulation  (4.11)  Is  obtained  by 

replacing  u and  v by  U and  V from  (4.13)  and  setting  W - S ■ ’J'.  . 

e e j 

This  technique  Is  simply  the  Galerkin  method  with  the  finite-element 


Interpolation  functions  as  trial  functions.  The  following  equations  are 
obtained  for  a single  element  0^: 

<»  - * r ^ - r ^ 

• <»  - ^ >ii“;  ^ ^ Sell”?;  * f ^ 

(B  + 2 2 °1J^  * 2 ®lj^l  ^ 2 ®ij^  ^ 2 ^ij^i 

■ <>  - + f - f 'ij”;  - f “ij*;  - f • 


where 
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“ -4 

e 

^ijk  ■ f{i^ 

°1J  ‘ 4^  ’''l.y'*'j 

®ij " fsi^  hjj,x  ‘**‘*y 

hi  - 4^  '"l.y’J'j.y  ‘»*‘»y  • 

We  note  that  In  early  phases  of  this  work  an  alternate  version  of  (4.11) 
was  used.  In  this  alternate  form,  the  first  term  on  the  right-hand  side 
In  (4.11)  was  Integrated  by  parts  to  produce  a nonconservative  form.  The 
nonconservative  finite-element  scheme  produced  an  unreasonable  answer. 

The  computed  velocity  profiles  were  In  error  at  and  downstream  of  the 
shock.  We  therefore  recommend  that  a weak  formulation  of  the  form  of 
(4.11)  be  used  rather  than  the  alternate  formulation  described  above. 

The  boundary  conditions  on  30  are  applied  In  the  standard  way  by 
setting  specified  model  values  of  velocity  on  the  boundary  either  to  be 
zero  or  to  take  on  prescribed  values  while  ellowlng  other  model  velocities 
on  the  boundary  to  be  unknowna.  This  method  works  satisfactorily  except 
at  the  leading  and  trailing  edge  of  the  airfoil.  Here  the  specified 
velocity  V Is  discontinuous.  Tor  example.  It  Is  zero  in  front  of  the 
airfoil  but  has  a nonzero  value  at  the  leading  edge.  The  finite  element 
model  naturally  assumes  a continuous  value  of  velocity  everywhere.  This 
dlscrepency  can  be  minimized  by  putting  In  very  small  elements  just 
before  the  leading  edge  and  just  after  the  trelllng  edge.  Then,  for 
example,  the  nodal  value  of  velocity  v on  the  leading  edge  In  the 
small  element  Is  given  a prescribed  value,  but  all  other  nodes  In  the 
small  element  and  on  3(2  are  set  to  zero.  Thus,  while  the  finite- 
element  velocity  model  s continuous.  It  approximates  s step  function. 
Hopefully,  the  small  elements  before  the  leading  edge  end  after  the 
trailing  edge  will  not  affect  the  conditioning  of  the  system. 

We  also  tried  an  alternate  procedure  In  this  work.  In  assembly 
of  the  elements  just  upstream  and  downstream  of  the  leading  and  trailing 
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edges,  a special  procedure  was  used.  The  compatibility  of  model  velocities 
V at  the  leading  edge  and  trailing  edge  was  never  enforced.  Then,  for 
example,  in  the  element  to  the  left  of  the  leading  edge,  the  nodal  value 
of  V at  the  leading  edge  was  aet  to  aero  while,  in  the  element  to  the 
right  of  the  leading  edge,  the  nodal  value  of  v at  the  leading  edge 
was  given  a prescribed  value.  This  method  did  not  produce  satisfactory 
results.  Apparently,  the  lack  of  compatibility  In  the  four  elements 
Involved  distorting  the  solution  near  the  leading  and  trailing  edges. 

Test  calculations  were  performed  to  determine  the  flow  about  a 
6-percent  parabolic  arc  airfoil.  The  mesh  used  Is  shown  In  figure  6. 

In  the  results  to  be  presented  here,  ■ 0.875  and  y ~ while 

. The  viscosity  parameter  a was  given  several 
values  to  demonstrate  the  convergence  of  the  solution  with  decreasing 
viscosity.  In  figure  3 the  approximate  distribution  on  the  airfoil 
Is  pictured  for  o 0.01.  The  finite-element  calculations  are  compared 
with  the  finite-difference  results  of  Murman  (1971,  1975)  for  K > 1.8  . 
Table  IV  shows  the  convergence  properties  of  the  solution  for  a ••  0.07  . 
The  polntwlse  value  of  model  velocity  u at  the  point  on  the  airfoil  of 
maximum  model  velocity  u of  the  converged  solution  Is  shown  for  the 
Iterative  steps.  A variable  time  step  was  used,  and  this  time  step  is 
Indicated  In  the  table.  Convergence  occurs  after  five  Iterative  steps. 
After  the  fifth  step,  the  algorithm  seems  to  oscillate  slightly  about 
the  solution.  Depending  on  how  strict  the  convergence  criterion  Is, 
this  oscillation  may  cause  an  Increased  number  of  Iterations.  Incor- 
poration of  a damping  term  Involving  a time  derivative  should  eliminate 
this  oscillation,  although  this  method  has  not  been  used  In  this  work. 


Figure  6.  Six  Percent  Perebollc  Arc  Airfoil 


TABLE  IV:  CONVERGENCE  PROPERTIES  OF  THE  APPROXIMATION  FOR  a » 0.07. 


Iteration 

Step 

u Velocity 

Component 

Time  Step 

1 

0.07714 

0.025 

2 

0.10730 

0.025 

3 

0.11309 

0.025 

4 

0.12995 

0.1 

5 

0.13476 

0.1 

6 

0.13409 

0.1 

7 

0.13461 

0.1 

9 

0.13481 

0.1 

11 

0.13496 

0.1 

Note: 


u velocity  components  represent  the  velocity  at  the  ninth  node 
from  the  leading  edge  on  the  airfoil. 


r 


4.6  Conclusions 

A hyperbolic  system  (In  time)  of  first-order  equations  Is  used 
to  calculate  both  steady  and  unsteady  flows  In  a manner  similar  to 
Magnus  and  Yoshlhara's  time-dependent  method.  First,  the  time  derivatives 
are  discretized  by  using  finite-differences  of  Implicit  Crank-Mlcholson 
type.  The  resulting  equations  are  linearized  with  Newton's  method. 
Regularization  techniques  are  needed  for  stability,  as  well  as  for  the 
treatment  of  shock  waves.  Different  artificial  viscosity  terms  are 
added.  Then,  the  space  derivatives  are  discretized  with  a Galerkln 
finite-element  procedure.  Stability  and  convergence  of  the  algorithm 
are  analyzed,  finite-element  Implementation  of  the  transonic  small- 
disturbance  case  Is  discussed,  and  a numerical  example  of  flow  around  a 
parabolic  arc  airfoil  Is  compared  to  Murman's  finite-difference  calcula- 
tions. Preliminary  results  Indicate  that  additional  work  Is  needed  to 
assess  the  validity,  applicability,  and  efficiency  of  the  algorithm. 
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APPENDIX  A;  PARAMETRIC  EXTRAPOLATION  METHODS  FOR  THE  SOLUTION  OF  THE 
HELMOHOLZ  AND  TRANSONIC  HARMONIC  PERTURBATION  EQUATIONS 


In  this  appendix  we  conelder  the  frequency  appearing  in  the 
perturbation  equation  as  a parameter  and  approximate  the  behavior  of  the 
solution  depending  on  this  parameter. 

Consider  the  one-dlmenslonal  exmaple: 


<t>^  + - 0 


(with  the  exact  solution  <J)  - sin  Uix) . 

For  small  values  of  U)  , the  solution  can  be  obtained  by  relaxation,  say 
and  If  we  assume  an  expansion  of  In  u , 

2 

♦ (u)  - + # (u  )(U)  - 0)  ) + ) 5-^  + .*.  , 


o'  - 2'  o''  2 

oco 


where 


3A  4>(wp  - 4)(ai3) 

((i)  ) S — - 

3(i)  ^ o'  Uj  - 


“o  - “2  ' 


“3  + “l 


(Wq)  =[<I'(u3)  - + ♦(Uj)]/[(a)j  - u^)/2]^  . 

Hence,  an  approximate  solution  4»  may  be  calculated  at  higher  values  of 
u)  In  terms  of  ()>(Uj)  , , and  ({>(0)3)  . 

A better  approximation  can  then  be  calculated  from  (*) 


- a)S  . 


For  this  example,  we  have 


*llotc  that  the  Iterative  procedure 


2.n 

♦xx  ♦ 


Is  divergent  for  u > u^j.  (since  It  can  be  described  by  the  time-dependent 
equation  + (»)^()>  - . 
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) 

I 

J 


2 I 

‘•'xx  “ ■ “ “o*  ■ * [co8(u)^x)](a)  - 0)^) 

+ x^  [sin  (oJjjX)]  (oj  - <*>q)^/2  I . 

Obviously,  the  polynomial  approximation  (Taylor  or  Chebyshev)  of  the 
sinusoidal  behavior  is  valid  with  some  restrictions,  which  may  not  ex- 
clude using  this  method  for  the  practical  application  of  relaxation  pro- 
grams beyond  the  critical  frequency. 

An  exponential  approximation  may  be  preferred;  namely  (Shanks  Formula): 

« 4>(Wq) 


(|)°  + 4)^  - 241^ 


, where  j4'  “ 4>(Wq)  + ~ ***0^ 


4)  « 4>(‘^q)  + ~ ^o^ 


1.  ....(w-w) 

+ 4>  (d)  ) o' 

o 2 

The  harmonic  (Fourier)  approximation  gives  the  exact  answer  In  this 
case  If  we  use  three  sets  of  data  to  determine  the  amplitude,  phase, 
and  frequency;  i.e.,  i^(aj)  “ A sin  (otojj^  + 0)  , 1 “ 1,2,3. 

For  the  two-dimensional,  transonic  harmonic  perturbation,  we  can  obtain  the 
solution  of  the  following  equation  for  different  (small  values  of  u by 
relaxation  techniques: 

(K  - 4)°)4)^  + 4)^  - (4)°  + 2lu))4)^  - 0 . 

' ^x  ^xx  yy  ^xx  x 

Then,  we  obtain  an  estimate  of  the  solution  4*  ^or  the  required  (high) 
value  of  b)  by  using  a polynomial,  exponential,  or  Fourier  approximation. 

The  final  answer  can  then  be  obtained  from  the  solution  of 


(K  - 4)°)4)^  + 4»^  - 4)°  4>^  - 210)4. 

X XX  yy  XX  x x 
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APPENDIX  B;  PROCEDURE  FOR  CALCULATING  PRESSURES  IN  THE  SHOCK  REGION 


Assume  form  of  solution 


♦ (x.y.t)  - (|>®(x,y)  + ee^^^^(x,y)  (A2.1) 

<l>jj(x»y.t)  - ♦“(x.y)  + ce*“^(|)^(x,y)  (A2.2) 

X°(y,t)  - X°(y)  + ee^4^(y)  (A2.3) 

Ue  assume  that  the  solutions  for  the  steady  problem  and  perturbation  problem 
are  knotm.  The  pressure  In  the  region  of  the  shock  must  now  be  calculated. 


Ue  consider  a point  A on  the  downstream  side  of  the  steady  shock.  If 
the  perturbed  shock  at  tine  t Is  upstream  of  X^  , the  velocity  Is  given 
as: 

- ♦°(X^,y)  + ce^*'<.^(X^,y)  (A2.4) 

for  X°  < X®  < X^  . 

If  the  perturbed  shock  at  time  t Is  downstream  of  X^  , the  velocity  Is 
given  as, 

♦jj(X^,y,t)  - ♦°(X°,y)  + ee^^'^^^CX^.y)  + (X^-X°)^^(X^,y)  (A2.5) 
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for  X°  < < X®  . 

Similarly,  for  a point  B on  the  upstream  side  of  the  steady  shock 

♦^(X®,y.t)  - 4.°(X®,y)  + ee^*=4.]^(X®,y)  (A2.6) 

for  X®  < X°  < X° 
and 

4.^(X®,y,t)  - ♦°(X°,y)  + ee^^<|.^(X°.y)  + (X®-X°)^^(x2,y)  (A2.7) 

for  X®  < X®  < X°  , 

where  X°  and  X^  denote  positions  on  the  upstream  side  of  the  shock  and 
the  downstream  side  of  the  shock,  respectively.  Thus,  to  calculate  the 
pressure  in  the  shock  region,  we  must  specify  a value  of  e so  that  we 
can  calculate  the  excursion  of  the  shock.  We  note  that  the  last  term  in 
eqs.  (A2.5)  and  (A2.7)  is  0(e). 

At  any  point  A (or  B)  the  pressure  Integrated  over  a cycle  is 
a weighted  sum  of  the  two  equations,  while  the  unsteady  pressure  is  given 
as: 

ce^‘(|)j[(X^,y)  for  X®  < X^  (A2.8a) 

♦°(X®,y)  - <^®(x^.y)  + ee^*'<>][(X°,y)  + (X^-X°)(>^(x^,y)  (A2.8b) 

for  3^  > X^  . 

Thus,  when  the  shock  passes  over  point  A,  there  is  a dramatic  jump  in 
the  pressure,  as  would  be  expected. 

Now,  we  compute  the  equivalent  unsteady  pressure  at  X^  by 
integrating  over  one  cycle  of  oscillation  of  period  T > ~ . 

T T T 

j ♦^(x^,y,t)dt  - /♦°(x^,y)dt  + e J e^%^(X^,y)dt 
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where  during  the  time  to  T2  the  shock  Is  downstream  of  point  A. 

Rewriting  this  equation,  we  get 


(|>®(X^,y)T  + e(f'(X^,y)  (4^)  - (|)°(X^,y)T  i 

♦ I ‘'‘A'y)  ♦ ( ^2  - ■'1)  i 

Note  that  j 

♦®(X^,y)  - ♦°(xj,y)  + , (A2.ll) 

so  the  middle  term  becomes  j 

I-  lo  ,0 1«2  - ’i*  • <"•'» 

Thus,  dividing  through  by  T , we  obtain  | 

! 

I 

♦;(A*.y)  - clf+ict^.y)  ■ ♦;(**.,)  - lA'l  X ( i 

o 
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and  we  see  that, 


- 12ir 


(^)ci 


(A2.13) 


♦“(x*.y)  - <a*.y)  - IKI 


(A2.14) 


- ‘I'^X^.y)  + [♦^(X^.y)  - <!>i(X^,y)  ] 


(A2.15) 


The  movement  of  the  shock  from  the  unsteady  perturbation  alters  the  mean 
or  steady-state  solution!  It  smooths  out  the  discontinuity.  A simple 
model  of  this  was  predicted  by  Tljdeman. 


steady  mean 


unsteady  shock 


In  addition,  the  unsteady  pressure  Is  altered  by  the  shock  motion. 

The  question  Is  whether  or  not  the  second  term  belongs  to 
the  unsteady  pressure  contribution.  It  Is  formally  of  order  (e°)  , 
but  It  arises  out  of  the  unsteady  motion.  We  conclude  that  the  lowest 
order  unsteady  solution  Is  of  order  (e°)  In  a narrow  region  surrounding 
the  steady-state  shock  position  and  of  thickness  0(e)  . 


APPENDIX  C;  HARMONIC  PERTURBATIONS  OF  UNSTEADY 
FULL  POTENTIAL  EQUATIONS 

In  this  appendix,  the  problem  of  an  oscillating  airfoil  Is 
analyzed  by  using  the  full  potential  equation. 


Governing  Equations 
A.  Quaslllnaar  Form: 


where 


(a^  - - 2uv^  + (a^  - 

XX  ’^xy  yy 

■ 2u(^  . + 2v^  ^ 

’^xt  ^yt  ^tt 

u ■ cos  a + 6 

X 


V - sin  a + * 

y 

^ 2 ^ [2(})^  + u^  + - l] 

M 

CO 

and  where  u,  v,  and  a have  bean  normalized  by  the  free-stream 
velocity  U^  . 

B.  Fully  Conservative  Form: 


■ If  “ I?  + If 

1 

where  p . [i  - + u^  + v^  - 1)]  ^ 

and  where  the  density  p has  been  normalized  by  the  free-steam 
density  p^^  . 
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Boundary  Conditions 

A.  F(x,  y,  t)  • 0 oi  y - B(x,  t)  - 0 


DF 


■ 0 or  uF  + vF  - 0 


Let 


Dt  X y 

(♦  )b  + sin  a - + cos  a]  • 


B(x,  t)  - B°(x)  + B^(x,  t)  , 


then  (♦  )g  + sin  a - [ ((J>^)g  + cos  o]  (B°  + B^)  + B^  ; 


moreover , 


<*y>B  " ^♦y^B® 


(<l>y)B  - 


B.  In  the  far  field,  the  mean  steady-state  value  of  (|)  Is  given  In 
terms  of  the  steady  circulation  r°  : 


i}»  - - tan  ^[\/l  - tan  (9  - o)j 


The  unsteady  far-fleld  formula  may  be  obtained  In  terms  of 
Klrchoff's  formula.  (Asymptotic  expansion  leads  to  a pulsating 
vortex  for  oscillation  airfoils.) 

Weak  Solutions 

The  jump  conditions  admitted  by  the  weak  solutions  are: 


A.  Shock  Relations 

B.  Contact  Surfaces  (Wakes) 


JpD^  + puD^  + pvD  I)  - 0 ; D(x,  y,  t)  - 0 


i 


(43.1) 


- l[pJI  - flpuj  + Dy  dpv]) 
D(x,  y;  t)  - X - x®(y;  t) 

|f-|lpD  - IpuJ  - |~IIpvJ 

(Conservation  of  mass  across  the  shock) 

From  Irrotatlonallty  conditions: 


F 


Hence 


: D : D - n(J)^ 
t X y u^t' 


jD  ! O'fryl 


(A3.1)  becomes 


where 


|5!  (Pfl  -lIPull  + 


,0 


is! 

3y 


cl 


(Tangential  velocity  along  the  shock  Is  preserved.) 

(B)  Across  the  wake,  pressure  Is  contlnous  and  so  Is  p ; l.e. 

p flDj.  + uD^  + vDyD  - 0 

(The  relative  velocity  normal  to  the  wake  surface  Is  zero.) 

OpI)  ■ 0- 

In  other  words,  |]p|J  - 0 or  [a  ])  - 0 

([24>^  + (<j)^  - cos  a)^  + (A  - sin  a)^|J  ■ 0 


til'*’  + 2 cos  a + (I<>yD 


+ 2 sin  a - 0 


Simplifications; 

1.  Assume  the  wake  vortlclcs  are  located  on  the  x axis. 

2.  Asssume  Is  contlnous  across  the  wake. 

The  Kutta  and  wake  boundary  conditions  thsJ  become 

O^t  D ■*•  B^x  ^♦x^^  - 0 ; X > 1 (y  - 0) 

(a  « 1)  . 
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An  Interesting  problem  would  be  to  trace  the  wake  and  to  use  both  relationships 
derived  before. 

Harmonic  Perturbations 

^ y)  + e Reje^^  <>^(x,  y)|  + • • • 

X - X°(x,  y)  + e Re|e^^  X' (y)  | + • • • 
where  the  boundary  condition  Is 

B(x,  t)  - B°(x)  + e Re|e^^  B^(x)  ] . 

If  we  subsltute  the  expansion  for  <l>  into  the  differential  equation  and  the 
jump  conditions  the  zero-order  problem  Is  exactly  the  steady  problem.  The 
first-order  problem  Is  the  perturbation,  where  X^  Is  the  shock  movement. 

A shock  fitting  method  for  Is  needed;  otherwise  we  may  use  a fully  conservative 
scheme  In  strained  coordlantes  as  In  the  small-disturbance  case. 


Zeroth-Order  Problem 

Differential  Equation; 
where 


(a°^  - u°^(^°  + 2u°v°<p°  + (a°^  - - 0 

^xx  '^xy  ' ^yy 

O . O 

u ■ cos  Qt  + A 

X 


u°  ■ sin  a + 4® 

y 


2 2 


-O^  _ 1 Y - 1 f o o 
a - ^ 1 u - V 


- 1 


Boundary  Conditions;  ® ®®*  “ ® 


Wake: 


•*yll  y - 0 “ ° » ll^®D  ■ ; X ^ 1 
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Far  Field: 


Jump  Conditions 


(J»  ■ - tan  ^ ( 'll  - tan  (0  - o)  J 

: DpV11.-(^)J[pV|1,  -» 


■11,-0 


Jameson's  Fully  Conservative  Schemes  or  rotated  schemes  with  shock 
fitting  may  be  used  to  solve  this  problem. 


First-Order  Problem 


Differential  Equation; 


(a°  - u°  )<t»^  + 2u°v'^(|)^  + (a°  - v°  )4)^  - N.H.  , 


where  N.H.  stands  for  nonhomegeneous  terms,  which  will  not  change  the  type  of 
the  equation. 


Boundary  Conditions: 


Wake: 


+ cos 


Far  Field: 


* , . o«>, . 0 ■ » 

■ 0 (for  the  present  time) 


Jump  Conditions:  1) 


x"  o 

X 


2)  !<.  |Ip°j/-DpVd  OpVD  OpVd  „ 

X X ^ I X X 


- loVl  - g- jpVl  . 

These  jump  conditions  must  be  Imposed  on  the  problem.  A shock  fitting 

procedure  is  necessary,  where  an  ordinary  differential  equation  in  x\  is 


*ln  this  case  multiply  the  equation  by  p°/(a°)^  . 


S8 


solved  simultaneously  with  the  partial  differential  equation  of  given 

and  X° 

An  alternative  to  the  approach  based  on  the  method  of  strained  coordinates 
may  be  used  here^  as  In  the  small-dlsturbance  case. 
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